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Abstract 

Using first principles molecular dynamics simulations, we have determined the threshold dis- 
placement energies and the associated created defects in cubic silicon carbide. Contrary to previous 
studies using classical molecular dynamics, we found values close to the experimental consensus, 
and also created defects in good agreement with recent works on interstitials stability in silicon 
carbide. We carefully investigated the limits of this approach. Our work shows that it is possible 
to calculate displacement energies with first principles accuracy in silicon carbide, and suggests 
that it may be also the case for other covalent materials. 
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Particle irradiation is a well known and extensively used technique, allowing to modify 
mechanical, magnetic, electrical and optical properties of materials. For instance, a suitable 
ion irradiation may harden a material, lead to a local oxydation state, or activate a magnetic 
order. The utility of ion irradiation is also well known for electronics, with the doping 
or gettering processes, and for radiation therapy. Besides, damage accumulation due to 
irradiation is also an important research field, related to space and nuclear applications. 

The interaction of an energetic ion with the matter is a complex phenomenon, especially 
at high energies. Impinging ions are simultaneously slowed down by inelastic collisions with 
electrons, and by elastic collisions with atoms. The displacement of lattice atoms leads to 
creation of defects and accumulation of damage. A key quantity, relevant to the process and 
different for each irradiated material, is the threshold displacement energy (E^). E^t may be 
defined as the minimal kinetic energy that has to be transferred to a lattice atom in order 
to create a stable Frenkel pair that survives at least 10^^^ s. For instance, E^ values are 
required as key input in large-scale irradiation simulation packages, such as SRIM/TRIM, 
extensively used for determining implantation profiles in doping processes, or for calculating 
damage accumulation in materials. 

This quantity is rather difficult to measure, since single created defects have to be identi- 
fied during experiments, and associated with a well-defined irradiation energy. Then, there 
has been an increasing number of works aiming at the E^ determination from molecular 
dynamics simulations. The procedure is simple: after a defined impulse given to an atom, 
which is usually called the primary knock-on atom (PKA), the evolution of the system is 
monitored. Once the transfered energy exceeds the Ed-, there is formation of a Frenkel pair 
in the system. As far as we know, all simulations but one Ij have been done with molecular 
dynamics and classical empirical potentials. In fact, several reasons hinder ab initio molec- 
ular dynamics. First, determining an energy threshold from molecular dynamics requires 
many runs, since the kinetic impulsion is progressively increased to find the threshold, and 
the procedure is stochastic due to non zero temperature. Second, large systems must be 
employed for high Ed. Finally, there is usually a Ed associated with each cristallographic 
direction and with each element in a multicomponent system, which considerably increases 
the number of runs. 

The relevant question is: are the empirical potentials precise enough to allow a correct 
determination of i?^? For metals, it seems that this is the case, embedded-atom potentials 
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being able to model metallic bonds with a good accuracy 

a a. 

There is less certainty 

for other materials, such as semiconductors or ceramics, since covalent bonds or charge 
transfers are hard to reproduce with potentials. In that case, an ab initio molecular dynamics 
determination would be extremelely useful. 

The silicon carbide is a good illustration of this issue. It is a promising material, with 
potential applications in electronics, as a replacement for silicon, and in nuclear technology. 
Silicon carbide is also very interesting from a fundamental point of view, since it can be 
considered as a model for zinc-blende two-component covalent materials. There have been 
several measurements of the Ed, with different techniques, but a large dispersion of values 
is obtained [4]. In lack of precise data, it is usually assumed that average values for C and 
Si sublattices are 20 eV and 35 eV, respectively. However, subsequent molecular dynamics 
studies did not clearly confirm these values. Average values were found from 17 to 40 eV 
"or C sublattice and from 42 to 57 eV for Si sublattice, with extreme values very different 
in add,t,o„, the „atu..e of the cheated defect. ,s d-ffe... f.cn one study 
to another. We have recently shown that these discrepancies are due to the use of different 
empirical potentials [l^. In fact, the kinetic energy required for the creation of a Frenkel pair 
is obviously related to the energy barrier that the lattice atom must overcome to reach an 
interstitial site. Empirical potentials usually give a poor description of these saddle states, 
especially for covalent materials. A precise ab initio determination would be invaluable in 
that case. 

In this paper, we report the first ab initio molecular dynamics determination of Ed- On 
the one hand, we show that such calculations are feasible, at least for covalent materials 
for which the vacancy-interstitial separation of the Frenkel pair is very small. On the other 
hand, Ed values have been obtained in /?-SiC for all high symmetry directions shown in 
figure [H for both Si and C lattices, with the first principles accuracy. Our results show that 
the use of available empirical potentials may lead to quantitative and qualitative errors, and 
that our calculated average values are close to the experimental consensus. 

The ab initio molecular dynamics calculations were performed using the plane-wave pseu- 
dopotential code GP based on the density functional theory (DFT) 12, 13]. The 



exchange-correlation potential proposed by Ceperley and Alder, and parametrized by Perdew 
and Zunger was used [l^. We considered a F-sampling of the Brillouin zone, and a 35 Ry 
kinetic cut-off. With those parameters, the calcultated lattice parameter oq = 4.34 A and 
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FIG. 1: Representation of the main crystallographic directions in /3-SiC. Carbon atoms are drawn 
in black, and silicon atoms in light grey (yellow in the electronic version). 

the bulk modulus 5 = 221 GPa were found to reproduce rather well experimental values, 
4.36 A and 224 GPa respectively [j. We also checked that pseudopotential cores did not 
overlap during simulations. All calculations were performed with a constant number of par- 
ticules, with a 64-atom cell (2ao x 2ao x 2ao), except for the Si PKA in the (100) direction 
where a 96-atom cell (3ao x 2ao x 2ao) was required to keep the PKA in the cell. A time 
step dt = 1 a.u was used during the ballistic phase of the simulation, then increased to 2 a.u 
during the relaxation phase. A thermostat was applied to recover the initial temperature 
of 300 K during the latter phase. The maximum duration of each run was 2.8 ps. If a 
stable Frenkel pair occured, the system was then completely relaxed to obtain the stable 
configuration. 

As an example, the figure [2] shows two possible cases in a typical threshold displacement 
energy determination, after a kinetic energy E is transferred to a silicon atom along the (111) 
direction. The PKA first moves from its equilibrium position along the (111) direction. If E 
is below the threshold displacement energy Ed-, in that case 22 eV, it returns to this location 
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FIG. 2: A Si PKA along the (111) direction. Carbon atoms are drawn in black, and silicon atoms 
in light grey (yellow in the electronic version). The silicon PKA is drawn in grey (orange in the 
electronic version), and the vacancy is represented by an open circle. A kinetic energy E is given 
to a Si atom, which is subsequently displaced. If E<£'rf, the PKA returns to its original location. If 
Fi>Eii, there is formation of a silicon vacancy Vsi and a silicon tetrahedral interstitial surrounded 
by four carbon atoms Sixc- 

and no Frenkel pair is created. On the contrary, if E is above Ed-, the PKA reaches an 
interstitial location in the lattice, leaving its original site free. Thus there will be formation 
of a Frenkel pair, i.e an interstitial and a vacancy, separated by a distance dpp- In this 
example, a vacancy and a silicon in a carbon tetrahedral site {Vsi+SiTc)) separated by a 
distance dpp = 0.87ao, are produced above the E^. 

There are several computational issues that are supposed to prevent the determination 
of Ed with first principles methods. Hence, the cell must be big enough to contain the 
PKA during all the simulation. Here, we have mainly used a 64-atom cell, which may be 
viewed as very small. However, in our simulations, the PKA does not move far away from 
its initial location before to be trapped in an interstitial site. Indeed in covalent materials, 
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and especially in ceramics, the vacancy-interstitial separation dpp is very short, often lower 
than Co- This is clearly in contrast to metals, for which dpp is several times qq. Also, the 
cell should be large enough to prevent cumbersome interactions between the PKA and the 
thermostat during the simulation. Hence, in silicon, it has been suggested that a 64-atom 
cell is too small with respect to this issue 16|. However we have recently shown that, in 
silicon carbide, the error due to the cell size problem is small compared to the discrepancy 



found between different calculation methods 



lOj. This is an important point, and we have 



performed an additional test with a larger cell (216-atom) and C(IOO) to check the validity 
of this assumption. We found no difference with the 64-atom cell, with a similar Ed value. 
Another issue is related to the time step. It must be small enough to insure the accuracy 
of atomic trajectories, especially during the ballistic phase of the simulation. Hence, we 
have used a time step of 1 a.u, so that the maximum displacement during one time step 
for a C PKA of 50 eV is less than 0.007 A, which is much lower than the upper threshold 
of 0.1 A recommended by Corrales et al. for low energy cascade events jl?]. Regarding all 
these points, we assert that the determination of by ab initio methods is feasible at least 
in ceramics, and, as it will be shown further, these calculations are required for determining 
accurately the threshold displacement energies and the created defects. 

We now describe and discuss our results. Table [J reports the calculated Ed values and 
the associated Frenkel pairs, obtained for PKA's on both C and Si sublattices in the main 
cristallographic directions. The corresponding defect configurations are reproduced in figure 
[3l Globally, our results show various dumbbells and Si interstitials in tetrahedral site Sire- 
For C[100] and an energy above 18 eV, the PKA recoils toward the nearest tetrahedral in- 
terstitial site and moves further until it forms a tilted CC(ioo) dumbbell interstitial with dpp 
equal to O.STcq. This configuration was previously described as the most stable CC dumb- 
bell |l^. Several CSi dumbbells were also identified. For C[110] and Ed equal to 14 eV, the 
C atom replaces its C first neighbor, which is subsequently displaced to create a CS'i^oio) 
with dpp = 0.47ao. This configuration is also found in the case of a Si PKA along the 
(110) direction, and an energy above 45 eV, with a different collision sequence. Considering 
now C[lll] direction, above 16 eV, the C atom heads for the tetrahedron defined by four Si 
atoms, and does not form a Ctsi tetrahedral interstitial as it could be primarly expected, 
but a slightly tilted CSi^oio) dumbbell with a Si atom. The Frenkel pair separation dpp is 
0.95ao- This is consistent with previous ab initio calculations from Lento et al., predicting 
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Direction 


(eV) 


Defect 


dpp («o) 


C[100] 


18 


Vn + tilted CC/^Qo) 


0.87 


C[110] 


14 


Vc + C5i(oio) 


0.48 


C[lll] 


38 


/ 


/ 


c[in] 


16 


Vc + CS'i (010) 


0.95 


C sublattice, 


weighted average: 19 eV 


Si[100] 


46 


Vsi + Sire 


1.52 


Si[110] 


45 


Vc + CS'i (010) 


0.48 


Si[lll] 


22 


Vsi + SlTC 


0.87 


Si[lll] 


21 


Vc + CS'i(oio) 


1.24 



Si sublattice, weighted average: 38 eV 



TABLE I: Threshold displacement energies in /3-SiC, calculated by DFT-LDA molecular dynam- 
ics, along the main crystallographic directions. The associated defects and resulting Frenkel pair 
separations dpp are also added. Vc, Vsi, CC, CSi and Sipc correspond respectively to a carbon 
vacancy, a silicon vacancy, a carbon-carbon dumbbell, a carbon-silicon dumbbell and a silicon in a 
carbon tetrahedral site. The average values are weighted for equivalent directions. For the C[lll] 
case, several defects were observed. 



the conversion of the Cpsi tetrahedral interstitial to the CS'i(oio) dumbbell interstitial jl9| . 
The last case for which a CSi dumbbell is obtained is the Si[III] with equal to 21 eV. 
Here the Si atom collides with its C first neighbor, displaces it, and returns to its original 
location. The resulting CS'i (oio) interstitial is separated from the vacancy by 1.24ao. Silicon 
tetrahedral interstitials surrounded by carbon atoms Sipci which were determined as the 
most stable tetrahedral interstitial 18,[l9|, [2^, were also created. The most simple case is 
Si [111] described in figure [2l Above 22 eV, the Si PKA directly moves toward the tetra- 
hedral site and forms a Sipc, 0.87ao away from the vacancy. A Si PKA along the (100) 
direction, with an energy higher than 46 eV, leads to the formation of a Sipc interstitial 
separated from the Si vacancy by 1.52ao, after a short collision sequence during which the 
Si PKA replaces another Si atom, this one moving in the following tetrahedral site. For 
the C[lll] case and an energy higher than 38 eV, several mechanisms, occuring for similar 
energies, were observed depending on the way the C PKA rebounded on its closest silicon 
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C[100] 



C[110] 




FIG. 3: Defect configurations for each considered crystallographic directions. Carbon atoms are 
drawn in black, and silicon atoms in light grey (yellow in the electronic version). Defects are 
drawn in grey (orange and red for C and Si atoms, respectively, in the electronic version), and the 
vacancies are represented by an open circle. 



neighbor. In the first mechanism, the C PKA rebounds without displacing the Si atom and 
forms CS'i(ioo), identical to C[110] and Si [110] cases. In the others, the C PKA encounters 
its Si first neighbor at short distance with enough energy to displace it to the next Sitc 
interstitial site. Afterwards, the C PKA sometimes returns to its original location, leading 
to a final configuration similar to the Si[lll] case, or it bounces backward, and after few re- 
combinations forms additional defects such as Csi antisite and carbon vacancy Vc-, as shown 
in figure O In this peculiar case, there is an uncertainty regarding the created defects, but 
for a similar Ed, a somewhat different result than in previous works 9|. Finally, regarding all 
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the different PKA's that have been studied, the created defects are always in fair agreement 
with the relative stability of defects found with static ab initio calculations. 

We have determined the average on both C and Si sublattices, by weighting each 
values of by the number of equivalent directions 2l|. Our average E^ are in very good 
agreement with the values usually considered by the fusion community: 19 eV against 20 eV 
for the C sublattice, and 38 eV against 35 eV for the Si sublattice. 

Several previous investigations were devoted to the determination of threshold displace- 
ment energies in silicon carbide with classical molecular dynamics, but with differences in 
calculated Ed values. Moreover, identified defects strongly diverged between all studies. 
Roughly it is possible to sortprevious results into two groups. The first one, related to the 
original Tersoff potential J, shows similar or slightly higher E^ values than in our work, 
but the created Frenkel pairs, mostly Vc + Crsi with low vacancy-interstitial separations, 
seems unphysical. It could be explained by the fact that the original Tersoff potential highly 



favored the formation of the Crsi interstitial [201,122]. Conversely, the second group, related 
to Tersoff potentials modified for short range interactions 6|, LZl, [sl, Isj , exhibits a more realistic 
defects production (essentially dumbbells), but also much higher Ed than in our work. Thus 
it seems hard to achieve a good description of the defect production together with accurate 
values of Ed using semi-empirical potentials. Ab initio tight-binding molecular dynamics 
have already been performed for two crystallographic directions [1], but the results were 
not convincing as the authors used a minimal basis set unable to take into account charge 
transfers. As a result, calculated defect formation energies, as well as Ed values, were not 
accurate. For example, the authors found Ed equal to 27.5 eV for a C PKA along the (100) 
direction, against 18 eV, here. State of the art first principles calculations are then required 
for determining accurately both Ed and formed defects. 

In conclusion, we have demonstrated that it is feasible to determine Ed in silicon carbide 
using ab initio molecular dynamics. This method could also be applied for silicon and other 
covalent materials. Average calculated values of Ed were found in very close agreement to 
the experimental consensus for both C and Si sublattices. Such an agreement, both with the 
experiment and the calculated defect formation energies, has never been found with semi- 
empirical potentials or tight-binding methods, and hence justifies the use of first principles 
methods for the determination of threshold displacement energies in covalent materials. 

This work was funded by the joint research program "ISMIR" between CEA and CNRS. 
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